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Abstract. We compute analytically the probability density function (pdf) of the 
largest eigenvalue A max in rotationally invariant Cauchy ensembles of N x N matrices. 
We consider unitary (/? = 2), orthogonal (J3 = 1) and symplectic ((3 = 4) ensembles 
of such heavy-tailed random matrices. We show that a central non-Gaussian regime 
for A max ~ O(N) is flanked by large deviation tails on both sides which we compute 
here exactly for any value of f3. By matching these tails with the central regime, we 
obtain the exact leading asymptotic behaviors of the pdf in the central regime, which 
generalizes the Tracy- Widom distribution known for Gaussian ensembles, both at small 
and large arguments and for any j3. Our analytical results are confirmed by numerical 
simulations. 
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Characterizing the distribution of extreme eigenvalues of random matrices is of 
paramount interest in many applications. A cornerstone in this field is the discovery of 
the Tracy- Widom (TW) laws [1] for the typical fluctuations of the largest eigenvalue of 
Gaussian random matrices. Let us first recall a few basic definitions and results for the 
Gaussian ensembles. Consider a N x N matrix X whose upper-triangular entries are 
drawn at random from a standard Gaussian distribution in the real (/3 = 1), complex 
(/3 = 2) or quaternion (/3 = 4) domain. Symmetrizing, we end up with a symmetric, 
Hermitian or quaternionic self-dual random matrix (respectively) whose N eigenvalues 
Ai, . . . , Ajv are real random variables, characterized by the joint probability distribution 
function (jpdf): 

iWAi, ■■■An)= B N (P)e-$^ A » J] \Xj - \ k f , (1) 

j<k 

where the normalization constant B^(P) is known from the celebrated Selberg's integral. 
This joint law allows one to interpret the Aj's as the positions of charged particles 
repelling each other via a 2d-Coulomb (logarithmic) potential: they are confined on a 
Id line and each one is subject to a harmonic potential. Here we focus on the largest 
eigenvalue A max = maxi<j<jvAj. What can we say about its statistical properties? This 
is a non-trivial question because, due to the all-to-all interaction term Ylj<k ~ ^kl 13 , 
the standard results about extreme value statistics of independent random variables [2] 
(namely, the existence of just three universality classes Gumbel, Frechet or Weibull) are 
no longer applicable. 

The average location of the largest eigenvalue is readily determined by the shape 
of the average density of the eigenvalues in the large N limit. For a Gaussian 
random matrix of large size N, the average density of eigenvalues (normalized to 
unity) Piv(A) = (^^j^(A — \)) has a semi-circular shape on the compact support 
[—\/2N, \/2N] called the Wigner semicircle 

Piv(A) « -j^ T Pw ( -7= J with p w {x) = -V2 - x 2 . (2) 

It thus follows that the average location of the largest eigenvalue is given for large N 
by the upper edge of the density support: 

(A max ) « V2N . (3) 

However, the largest eigenvalue fluctuates from one realization of the matrix to another. 
What can be said about the full probability density of A max ? 

From the jpdf Q it is easy to write the cumulative distribution of A max as a multiple 
integral: 

F(A max <w) = B N {(3) Yl / d\ F jato (Ai, ...,\ N ), (4) 

i=i 

which can be interpreted as the partition function of a Coulomb gas in the presence of 
a hard wall in w. Carrying out this multiple integration is a non-trivial task. It turns 
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out that it is necessary to deal separately with two natural scales for the fluctuations 
in the asymptotic limit. Typical (small) fluctuations scale as ~ iV -1 / 6 , while atypical 
(large) fluctuations scale as iV 1//2 , and the two corresponding distributions are described 
by different functional forms: 

• Typical fluctuations: Forrester j3] followed by Tracy and Widom [dj realized that 
in the large N limit, the largest eigenvalue follows the law 



/l max ~ 



(5) 



with a\fl = l/\/2 and a± = 2 _7//6 and where the random variable xp has an N- 
independent distribution, F(xp < x) — J~p{x) called the Tracy- Widom distribution 
PQ, which has highly asymmetric tails: 



T'Jx) ~ exp 



exp 



24' 

2(3 
— —x 



3/2 



as x — > —00 



as x — > 00 . 



(6) 
(7) 



These TW distributions also describe the top eigenvalue statistics of large real 
[U [5] and complex [6] Gaussian Wishart matrices, which play an important role 
in Principal Component Analysis of large datasets. Amazingly, the same TW 
distributions have emerged in a number of a priori unrelated problems [7J such as 
the longest increasing subsequence of random permutations [S], directed polymers 
[6j [9] and growth models [10] in the Kardar-Parisi-Zhang universality class in 1+1 
dimensions, sequence alignment problems [UJ, mesoscopic fluctuations in quantum 
dots [12J, height fluctuations of non-intersecting Brownian motions [321 [Hj and 
also in finance [15]. Remarkably, the TW distributions associated to the Gaussian 
Unitary and Orthogonal Ensembles have been recently observed in experiments on 
nematic liquid crystals [16J. 

Atypical fluctuations: the TW laws do not account for atypically large 
fluctuations, e.g. of order 0(y/N) around the mean value y2N. Questions related 
to large deviations of extreme eigenvalues have recently emerged in cosmology [17] 
and disordered systems [TBI US]) an d in the assessment of the efficiency of data 
compression |20j . Recently, the large deviations of the largest eigenvalue of Wishart 
matrices have been measured in experiments involving coupled fiber lasers |21j. The 
probability of atypical large fluctuations, to leading order for large N, is usually 
described by two large deviation (or rate) functions ip-(x) (for fluctuations to the 
left of the mean) and ip+(x) (for fluctuations to the right of the mean), in such a 
way that the probability density of the largest eigenvalue reads: 



d 

dw 



exp 



P(A max < w) « { 



1 XV / v 



VN 
2N 



exp 



w < 



w > 



^2N&:\w- V2N\ 
^2Nk\w - V2N\ 



t 0(y/N) 

0(JV-V6) 

(8) 
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Note that while the TW distribution J-p(x), describing the central part of the 
probability distribution of A max , depends explicitly on (3, the two leading order 
rate functions ip T (z) are independent of (3. Exploiting a simple physical method 
based on a Coulomb gas analogy (see below), the left rate function ip-(z) was first 
explicitly computed in [TBI US], while the right rate function ip + (z) was computed 
in [20]. A more complicated, albeit mathematically rigorous derivation (but only 
valid for /3 = 1) of ip+(z) in the context of spin glass models can be found in [22]- 
More recently, the various subleading corrections to the leading behavior have been 
explicitly computed using more sophisticated methods both for the left tail [23j, as 
well as for the right tail [2H . Note also that large deviations for the smallest 
eigenvalue were also studied for Wishart [26J and Jacobi [27] ensembles. 
The physical mechanism responsible for the left tail is very different from the one 
on the right. Computing the probability distribution P(A max < w) of the top 
eigenvalue is equivalent to computing the free energy of an interacting Coulomb 
gas in presence of a hard wall at w [HI [19]. For the left tail of P(A max < w), with 
w < \/2N, the charge density is pushed by the wall, which leads to a complete 
reorganization of all the N charges and thus costs an energy difference, compared 
to the Wigner sea (|2]), of order 0(N 2 ). In contrast, for the right tail, the dominant 
fluctuations are caused by pulling a single charge away from this Wigner sea and 
the energy difference is only of order O(N) [20]. The central TW region of the 
distribution in Eq. Q, describing typical fluctuations, matches smoothly with the 
two tail behaviors at its flanks [HI [19], EQ] . 

Given the apparent robustness of the predictions stemming from random matrices 
with Gaussian independent entries, it is of paramount interest to investigate whether 
these TW distributions and their large deviation tails for A max actually hold for a 
broader class of matrices. Results are available for 1) non-invariant ensembles, with 
independent and identically distributed (i.i.d.) entries, and 2) invariant ensembles, and 
are summarized here: 

1) It is known that TW holds asymptotically for symmetric N x N matrices with 
i.i.d. entries of variance 1/N, such that all moments are finite [28J. On the other 
hand, when the distribution of i.i.d. entries decays as a power law, ~ |Mjj| _1_At , 
the case \i — 4 leads to a new class of limiting distribution, while for fi > 4 the 
TW still holds asymptotically and for ji < 4 the statistics of the largest eigenvalue 
is governed by a Frechet law [151 [2S1 ES ED] • 

2) Much fewer results are known for rotationally invariant ensembles. The TW law has 
been established for the classical Wishart ensembles in [H [51 E] and for the Jacobi 
ensemble in [32J. Disordered ensembles of random matrices were studied in |33j 
where continuous transitions between TW and other extreme value distributions 
were found. In the context of Levy stable ensembles, the largest eigenvalue 
distribution for the so-called Levy-Smirnov ensemble has been derived in [34] . 

Here we address these challenging questions on top eigenvalue statistics beyond TW 
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and focus on yet another instance of random matrices which incorporates (i) highly non- 
Gaussian statistics and (ii) correlated entries, while retaining rotational invariance. We 
consider the Cauchy ensemble oi N x N matrix H, which might be symmetric {(3 = 1), 
Hermitian {(3 = 2) or quaternionic (f3 = 4). The weight associated to H is given by 

P(H)oc [det(l* + H 2 )]-*"- 1)/2 - 1 , (9) 

where 1 at denotes the N x N identity matrix. Note that Eq. ^ is manifestly invariant 
under the similarity transformation H — > UHU \ with U an orthogonal {(3 = 1), 
unitary {(3 = 2) or symplectic (/3 = 4) matrix. Due to this invariance, the jpdf of the N 
real eigenvalues can be straightforwardly written as: 



J°int(Al, • • • , X N ) OC Yl n A 2^(JV-l)/2+l II \ X 3 ~ X kf ■ ( 10 ) 

i=l 1 i> j<k 



As in the Gaussian case, this expression (10) allows to interpret the Aj's as the positions 



of charged particles (with say positive unit charge) repelling each other via the 2d- 
Coulomb (logarithmic) interaction. Here they are confined on the real line and interact 
in addition with a fixed particle, with charge — (N — 1 + 2//?), which is placed at the 
point of coordinate (0, 1) in the complex plane. 

This ensemble has been studied in the literature in Ref. [35J in the context of 
mesoscopic transport as a model of a quantum dot which is coupled to the outside world 
by non ideal leads containing N scattering channels. It is also one of the paradigmatic 
examples of application of free probability theory in the context of random matrices 
[36| [37]. A remarkable difference with the standard Gaussian ensemble is that the 
density pjv(A) is independent of N, Pat(A) = p*(A) (for all values of N), and its support 
is unbounded and given, for any (3, by [381 EH] 

^ A ) = -7X^' AG (-00,00). (11) 

7T 1 + A z 

In this paper, we study the statistics of the largest eigenvalue A max for the Cauchy 



ensemble in Eq. (10) for any (3, while up to now only the case (3 = 2 was considered 
in Ref. [101 SI] - and even in this case the large deviations of A max have not been 
computed. To estimate the typical scale of A max , denoted by A max , we note that for a 
general matrix model it satisfies 

POO 

/ p N (X)d\ w 1/N , (12) 
as the fraction of eigenvalues to the right of the maximum (including itself) is typically 



1/N. Substituting pjv(A) = p*(A), from (11), in (12), one obtains that A max ~ O(N). 



As we show below, this gives rise to three distinct regimes in the fluctuations of A r 
[as in Eq. (pi) for Gaussian ensembles, where in that case A max ~ 0(vN)]. 
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Summary of new results: Before presenting the details of the calculations, it is 
useful first to summarize our main results. We show that the pdf of A max , P(w,N) = 
<9u,P(A max < w), displays three different regimes to leading order for large N 

exp[-(3N 2 tfj(w)] w < iV 
P(w,N)^\ N- X fp{w/N) w~N (13) 
N(j)(w) w^> N 



log(2) 



where ip(w) and <fi(w) are independent of (5 and are given by 

V>H = \u lo § K + l ) ~ lo § { w + v^TT)) + ^f 1 , (14) 

and 

<f>(w) = — 2 . (15) 

7TW 

The full expression of the scaling function fp(x), describing typical fluctuations and 
playing the role analogous to the TW distributions in the Gaussian case ([8]), can be 
computed explicitly only for j3 = 2 J^Qj, E] , where it can be expressed in terms of the 



solution of a Painleve-V equation [see Eqs (48), (49) below]. The function f^ix) has its 



support over [0, oo). From the Painleve-V equation, we can check explicitly that f2(x) 
has the following asymptotic tails: ^(z) ~ exp (— l/8x 2 ) as x — » and f2^x) ~ 1/nx 2 
as x — > oo. While we were unable to compute fp(x) analytically for arbitrary 0, we are 
however able to predict its precise tails for any (3 > 

{exp f-T7-j) , x -»■ 
i V 16a; 2 7 (16) 
— - , X -> oo . 

We achieve this by smoothly matching the central regime A max ~ O(N) with the left 
(Amax *C iV) and the right (A max ^> iV) large deviation tails that we can compute for 
arbitrary 0. 

By comparing the results for the Gaussian case in Eq. (JsT) and the Cauchy case 



in Eq. (13), we see that the behavior for the right tail are rather different in the two 
cases. Indeed, quite generally, one can show (see later) that far to the right of the 
central peak, the pdf P(w, N) can be simply expressed in terms of the average density 
of states as P(w,N) TVpjv(A). In the Gaussian case, Pat(A) has a bounded support 
in the N — > oo limit (Wigner semi-circle). However, for finite but large N, there is a 
nontrivial correction to the density arising from the vicinity of the edge of the semi- 
circle [13], which leads to the nontrivial right large deviation behavior described by 
ip+(z) in Eq. ([§]). In contrast, in the Cauchy case, the average density of states has an 
unbounded support and there does not seem to be any nontrivial edge corrections for 
finite iV as in the Gaussian case. 



Right large deviation tail: Let us start by deriving the expression of the pdf P(w, N) 
for the right tail, i.e. where w 3> N. For this purpose, we define, for each eigenvalue 
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Aj, a binary variable <7j such that Cj = 1 if A, > w and o~i = if Aj < u>. Then, the 
probability that the region [w, oo) is free of eigenvalues can be written as 

P(A max <w) = ([1 - ai)[l -a 2 ) ■■■[!- a N )) , (17) 



where the average (•) is over the jpdf of the eigenvalues (10). Expanding the product in 



(17), one gets 

POO 

P(A max < w) = 1 — N / pjv(A)dA + two-point + three-point + . . . , (18) 

J w 

where 'two-point' means a double integral involving two-point correlation function (and 
similarly for 'three-point' etc). When w — > oo, keeping N fixed, which corresponds 
to the extreme right tail, one can show that all higher order contributions vanish 



Hence, substituting pjv(A) = p*(A), from (11), in (18) one obtains 



oo 



F(A max < u;) ^ 1 - iV / p*(X)d\. (19) 



w 



Taking the derivative of Eq. (19) with respect to (w.r.t.) w, and using p*(A) ~ l/(7rA ) 



for A 3> 1 (11), yield the result announced in Eq. (13) for w ^> N and <p(w 



l/(nw 2 ) (15). 



Left large deviation tail: Let us now focus on the left tail w N, which we tackle 
using a Coulomb gas technique. By definition: 

pia„„ < w] = j^'S'^M! 1 ?! ■ < 2 °> 

^(-oo.oo)^ 1 li=l dAi nomt[M, • • • , X N ) 

We can represent Pj mt(Ai, . . . , Ajv) in the Boltzmann form, Pj mt(Ai, . . . , Ajv) oc 
exp(— (j3/2)E[{X}]) where the energy function is given by: 

N 

Ail- (21) 



E[{\}] = ( N -l + ?)± HI + A?) - J> |A 



This thermodynamical analogy, originally due to Dyson [44j, allows to treat the system 
of eigenvalues as a 2d gas of charged particles (Coulomb gas) confined on the real line, in 
equilibrium under competing interactions. In the large iV limit, the Coulomb gas with 
iV discrete charges becomes a continuous gas which can be described by a continuum 
(normalized to unity) density function g w (x) = (^/N)J2i^( x ~ \-). Consequently, 



one can replace the original multiple integral in (20) by a functional integral over the 
space of Q w {x). This procedure, originally introduced by Dyson [B] for the problem 
without a wall, was first used successfully for the Gaussian case with a wall in [18], [19] 
and subsequently was found useful in a number of different contexts (see for example 
[45| HQ H7] and references therein). Introducing the constrained density of eigenvalues 
g w {x), one has 

f ( I3N 2 \ 
<w]oc V[g w ]expl — S w [g w ] J , (22) 



max 
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where the action is given by: 



ln(l + x 2 )g w (x)dx — 



— oo J — oo 



c 



1 



dx dx' q w (x)q w (x') In \x — x'\ 

(23) 



where C is a Lagrange multiplier enforcing the normalization of the density to 1. 
The cumulative distribution in Eq. (22) must be normalized, i.e., the proportionality 
constant in Eq. (22) is determined from the condition: P[A max < w] — > 1 as w — > oo. 
The physical meaning of the action <Su,[£> TO ] is readily understood as the free energy of 
the Coulomb gas, whose particles are constrained to lie to the left of a hard wall at 
x = w. The fact that the free energy is proportional to iV 2 (and not just N) is a 
consequence of the strong all-to-all interactions among the particles: the free energy 
of this correlated gas is dominated by the energetic component ~ 0(N 2 ), while the 
entropic term [~ O(N)} is subdominant in the large N limit (recently this entropic 
term has been computed explicitly for the Gaussian [IB] and the Wishart-Laguerre 
ensemble [IS]). For large N, the functional integral in Eq. (22) can be evaluated by the 



saddle point method. Clearly, the equilibrium configuration of the gas (the saddle point 
density) in the presence of the hard wall at x = w will be determined by minimizing 
the action in (23). Following the result in the Gaussian case [HI EES], we can anticipate 



that the modified equilibrium (saddle point) density q* w (x) in presence of the wall will 
display an integrable divergence as x — > w~ . In terms of the equilibrium density o^(x), 



the probability of the largest eigenvalue reads from (22): 



ax < w] « exp [-/3A^H] , 1>(w) = - (S w [q* w ] - Soof&J) . (24) 
The additional term ^oof^j comes from the normalization constant. In order to 



compute £?* (x), we vary the action 



<V5 



0, yielding: 



/w 
dx'g 
-oo 



' „* t~.l 

w 



x) In \x — x'\ , x G (— oo, w] . 



(25) 



Differentiating (25) w.r.t. x and setting x' = w — r and x = w — z yields: 

w — z 



Pr 



dr 



Qw(r) 



(26) 



/ T — Z 1 + [W — Z) 2 

where Pr stands for Cauchy's principal part and the shifted density q w {t) = g^(w — r 



is introduced. The inversion of the half-Hilbert transform (26) is the main technical 



challenge. However, the task can be fully accomplished [501 EI], and the general solution 



of the equation (26) of the Tricomi type can be written as: 

. , > 1 



Qw{T 



TT 2 \/T 



Pr 



ds 



1 y/s(lV - S) 

S — T 1 + (w — s)'' 



+ B )■ , r > 0, 



(27) 
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Figure 1. Plot of g^(x) = g w (w 

compared to the results of our numerical simulations (symbols 
N = 100. 



x) in Eq. (28) for different values of w (solid lines), 

for which (3 = 2 and 



where B is a constant enforcing the normalization J °° drg w (r) 
straightforward manipulations one finds that B = and eventually 



Qw(r) 



(k w +w) 1 / 2 -( W -T)(k w -w) 1 / 2 

1+(w-t) 2 



A^Z ) ((T^-fc^)(^-r)+fc tu (w-r) 2 +fc w )(fc m +t») 1 / 2 +r(t»-r)(fc tu --w) 1 / 2 
" X ~' ' k w [(w-rf+l] 



w>0 
w < 



After 



(28) 



where k w = \/w 2 + 1. One can check that J °° drg w (r) = 1 and g w (w — r) ~ ^n^2\ f° r 
w — > oo, which yields the result announced above ( JTlj ). In Fig. [TJ we show a plot of 
£>* (x) for different values of w, and compare our exact analytical expression to numerical 
simulations (see below for more details about them). 

Evaluation of the saddle point action: From the solution f?*(x) (depending 



parametrically on w) we can evaluate the action (23) as: 



/w pw r-w 

ln(l + x 2 )g* w (x)dx - / / dx dx' g* w {x)q* w (x') In \x - x 
oo J — oo J — oo 



(29) 



The double integral in (|29|) can be written in terms of a single integral, using the 

(30) 



following trick. We multiply (25) by g* w (x) and integrate over x to get eventually: 

C 
2~ ' 



/w riv ^ f>w 

I dx dx'g k w {x)g k w (x r ) In \x - x'\ = - / dxg* w (x) ln(l + x 2 ) + 
-OO t/ — oo J — oo 



' —CO J — oo 

Therefore, the action reads: 



1 /' u ' 



2~ 



(31) 



Given the different expression of g w (r) for w > and w < (28), these two cases have 
to be treated separately. We first consider the case w > where the constant C can be 



determined by setting x = in (|25|), yielding: 

C = 2 I dxlia \x\g* w (x) 



(32) 
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(33) 



Let us focus on the action «S W [^>*] in (33). Setting x = w — r and recalling that 



Q w {t) = Qw(w — r), we can rewrite the action, for w > 0, as: 

roc /*oo 

S w [q w ] = 7; \n (l + (w - t) 2 ) g w (r)dr - dr\n 
* Jo Jo 



W - T\Q W [T) , 



(34) 



where g w (r) is given in Eq. (28). After cumbersome manipulations, these integrals can 
be computed explicitly and we finally obtain a remarkably simple expression: 



S w [g w ) = i Q log (w 2 + 1) - log (w + Vw 2 + 1^ + 



31og(2) 



In particular, it admits the large w expansion 

1 



log 2 + 



+ 0(w~ 



(35) 



(36) 



For w < 0, we start again with the expression in Eq. (31) but in this case C can 



not be evaluated by setting x = in Eq. (25) as x = is not in the support on (x 



see Eq. (25)]. One can however evaluate C by setting x = w in Eq. (25) to get 



C 



Q* w (x) log \w — x\dx — log (1 + w z 



(37) 



After some manipulations, one can finally show that, for w < 0, the expression of S w [g n 



is also given by Eq. (35). Therefore, as announced in Eq. (13), the right tail of the 



cumulative distribution of A max is given by 
« < w] « exp (-/3iV 2 </>(uO) 

= ^ - <Soo[^]) = ig log (w 2 + l) - log (to + a/w 2 + 1^ 



+ 



2 

log(2) 



The rate function ip{w) has the following asymptotic tails 
ij){w) ~ l/[16iy 2 ] as»4 +oo , 



ln(|w|/2) + 0(\w 



-2^ 



as w — >■ — oo . 



(38) 

(39) 
(40) 



Substituting the asymptotic tail of ip(w) as w — > — oo in P[A max < w] ~ exp (—f3N 2 i(j(w)) 
then predicts a power law tail of the distribution as w — > — oo: P[A max < w] ~ 
[2/| W |]^ 2 / 2 . 
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For f3 — 2, it is interesting to compare the result (38 ) with previous works on the pdf 
nax , P(w,N) (HUH!]. In [10], the authors derived a nonlinear differential equation 



of A n 
for the function 

a(w,N) 



:i w 



P(w,N) 



[l + w 1 







lnP[A r 



< w] 



(41) 



f™ oo P(z,N)dz ^ ' ~ J dw 

They found that, for arbitrary positive N, a(w, N) = a satisfies a nonlinear differential 
equation 

(1 + w 2 ) 2 (a") 2 + 4(1 + w 2 ){a'f - 8wa(a') 2 + Aa 2 a' + 4N 2 (a') 2 = , (42) 

with the asymptotic property 

a(w, N) , w ->• oo . (43) 

Using the results of Ref. [12], it was shown in Ref. [ID] that Eq. (42) can be transformed 
into a Painleve-VI equation. 

Let us now check if our result in Eq. (38), after setting (3 = 2, is compatible with 
the above differential equation in the large N limit. Let us first evaluate a(w, N) defined 
in Eq. (41) by using our large N prediction in Eq. (38). Evidently, for large N but 

(44) 



finite w, we anticipate that a{w) must be proportional to N 2 , i.e., 

a(w,N) -> N 2 cr(w) , 
where the function a{w) is independent of iV and our result predicts that 

-(i + w 2 )i)\w) 



-(vT+ 



(45) 



On the other hand, substituting a(w,N) = N a(w) in the differential equation (42) 
and keeping only the leading order 0(N 6 ) terms, one finds that a(w) must satisfy the 
following nonlinear differential equation 

(1 + w 2 ) (a 1 ) 2 - 2waa' + a 2 + a' = . (46) 

Hence, in order to be compatible, our predicted form a(w) = |(a/1 + w 2 — w) must 
satisfy this nonlinear differential equation (46). Indeed, one can check quite easily that 
&(w) = |(\/l + w 2 — w) does satisfy Eq. (46), with the correct boundary condition (43): 
this thus provides a non-trivial check of our result for the rate function ip{w) (|38j). 

The central part: We now come to the central part of the pdf of A max where 
Amax ~ 0(N). From standard scaling arguments one expects, for w ~ O(N) 



(47) 



for any value of (3, as announced in (13). For (3 = 2, one can check that it is consistent 
with the large iV limit of (42) from which one obtains 

d 



hi?) 



dx 



exp 



1/x , .dz 
r{z)- 



x > 



(4£ 
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where t{x) = r satisfies a Painleve V equation on the interval [0, +00) 

x \ T y + 4(r - xr') (r - r\x + r')) = . 



(49) 



Interestingly, the same equation (|49j) also appears in the expression of the pdf of the 
(scaled) spacing between consecutive eigenvalues in the bulk of the spectrum of GUE 
random matrices [52]. Hence one expects that f 2 {x) can also be expressed in terms of 
a Fredholm determinant involving a kernel very similar to the sine-kernel [2HJ E3] • This 



can be shown directly from Eqs. (10) and (20), using orthogonal polynomials together 



with the asymptotic analysis performed in Ref. [31] to obtain [4"T] 

1 sin( l/yi - I/1/2 ) 

7T 



f 2 {x) = d x det(I - U X KU X ) , K( yi ,y 2 ) 



(50) 



yi - V2 

where I denotes the identity operator and the operator EL,, is the projector on the interval 



[x, +00). In Eq. (50), it is easy to see that under the change of variable y — > l/(iry) } 
K(y 1 ,y 2 ) translates into the standard sine-kernel [331 139"] . 

What are the asymptotic behaviors of f 2 (x) ? From Eq. (49), it is easy to see that 
t(x) ~ — x 2 /A, as x — > 00. This yields the small argument behavior of f 2 (x), from (48): 

1 

~8x? 



\nf 2 (x) 



x -»■ . 



(51) 



On the other hand, from Eq. (49), one obtains the small argument behavior of 



t(x) ~ —Xx, when x — >■ 0, but the constant A remains unspecified by this equation. 
Its determination can however be obtained by analyzing the Fredholm determinant in 
Eq. (50) which yields t(x) ~ —x/n, for x — > [32], and finally 

1 



h{x) 



x — > OO 



7rar 



(52) 



As we show now, these asymptotic behaviors for j3 = 2 in Eqs. (51) and (52) 



can be obtained straightforwardly, and generalized to any value of f3, by exploiting the 
matching between the central part and the right and left tails of the distribution of 
A max . Indeed, if we study the distribution of A max when A max approaches its typical 



value A max ~ 0(A) from below, 
above (38) that, given that ip{w) • 

P(A max < w) ~ exp 



i.e. from the left tail, one expects from our result 
- l/(16w 2 ) for w > 1 



-N 



P 



16w 2 



,w > 1 & w < N . 



(53) 



This expression (J53j) is obviously a function of the scaled variable x = w/N and coincides 
with the small argument behavior of the pdf fp(x) when x — w/N 1. This yields the 



asymptotic behavior announced in the introduction (16) and coincides, in particular, 
with the above result for j3 = 2 (51). 

Similarly, if one studies the pdf of A max when A max approaches its typical value 
A max ~ O(N) from above, i.e. from the right tail, one expects from our result in (13) 
that the pdf behaves like 

N 



P(w,N) 



7r«r 



w > N , 



(54) 
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co/N 



Figure 2. The scaling relation P(w,N) = N^ 1 f 2 (w /N) for /3 = 2 is tested against 
numerical simulations varying N. Montecarlo data are compared with the solution of 
Eq. (49) showing perfect agreement. Inset: the asymptotic behavior of the tails (16) 
is verified. 



which is obviously a function of the scaled variable w/N times 1/N. If one assumes 
that this behavior matches with the large argument behavior of the central part where 
P(w,N) ~ (1/N) f 13(10/ N) where w 3> N, one obtains that fp(x) ~ 1/(tcx 2 ) for large x, 
as announced in (fl3l). 



Numerical simulations: We now present the results of our numerical simulations. 



The distribution of the eigenvalues (10) can be directly simulated by exploiting the 



thermodynamical analogy, interpreting the eigenvalues Aj as a one- dimensional gas 
of charged particles with Coulomb interactions [3H1 Hence the distribution of 

the maximum eigenvalue can be obtained by means of a Metropolis algorithm, the 
acceptance rate between two different configurations being e~^ AE where the energy 



function i?[{A}] is given by Eq. (21). Within this numerical scheme, we have tested 



our analytical predictions in Eqs. (13), (16). In Fig. [2] we show a plot of NP(w,N) 
as a function of the scaled variable x = w/N for /3 = 2 and different values of 
iV = 80, 160, 320. The good collapse of the different curves is in agreement with the 



scaling form predicted for the typical fluctuations in the central regime (13), (47). The 



solid line corresponds to a numerical evaluation of fi(x) in Eq. (48) obtained by solving 



numerically the equation for r(x) in Eq. (49) (together with the asymptotic behavior 
t(x) ~ —x/tt for x — > 0): the agreement with the numerical data generated by the 
Metropolis algorithm is excellent. In the inset of Fig. [2j we see that our exact results 



in Eq. (16) describe very well the asymptotic behaviors of f2(x). In addition, in Fig. [3J 
we show a plot of — In P(w, N) as a function of w for iV fixed and large (N = 100) and 
for different values of f3 = 1, 2 and 4. This plot demonstrates the existence of the three 
distinct regimes as predicted in Eq. (Tl3|) . It shows a very good quantitative agreement 



with the large deviation functions ip(w) (38) and 4>(w) (15) which we predicted above. 
Note that ijj(w) has a support on the full real line, not only on the positive axis. However, 
for w < 0, P(w, N) is extremely small for N = 100 and is thus extremely hard to 
compute numerically (see Fig. [3]). Simulations for smaller values of N, not shown 
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Figure 3. Montccarlo results for the Coulomb gas system. The logarithm of the 
probability density is shown in the case N = 100 for j3 — 1 in a), (3 = 2 in b) and 
(3 = 4 in c). We compare it with our predictions for the left ( [38] ) and right (15) tails 
in dashed line. For (3 — 2 we also show the results for f2(x) ( |48[ ), in solid line, which 
describes the typical fluctuations. Note that one expects that tp(w) has a support on 
the full real axis, including w < 0, where P(w, N) is however extremely small and thus 
very hard to measure numerically for N = 100. 



here, agree with our calculations for w < 0, although they suffer from strong finite N 
corrections. 

In summary, we have investigated the fluctuations of the largest eigenvalue A max 
oi N x N random matrices belonging to rotationally invariant Cauchy ensembles. We 
have considered symmetric (/3 = 1), Hermitian (/3 = 2) and quaternionic (/3 = 4) 
matrix ensembles. We have identified three different regimes for the pdf of A max : (i) 
Amax "C N (left tail), (ii) A max ~ O(N) (central part) and (iii) A max 3> N (right tail). 
We have obtained exact results for the large deviation tails, both left and right and for 
any j3, which allows to obtain also the leading asymptotic behaviors of the pdf in the 
central regime, which generalizes the TW distributions known for Gaussian ensembles. 
These exact results have been confirmed by thorough numerical simulations. The exact 
expression for the central part of the distribution, describing the typical fluctuations of 
Amax ~ 0(N), beyond the case (3 = 2, remains a challenging problem. 
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